perm filename MFRAST.SAI[MF,DEK]17 blob sn#731855 filedate 1983-11-26 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00011 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	entry begin comment The semantics module of METAFONT.
C00006 00003	Procedures for generating and storing pens
C00041 00004	Procedures for plotting points on the raster
C00057 00005	Beginning of the main procedure \\{drawit}
C00065 00006	The routine that plots a cubic
C00079 00007	procedure ccubics # compute the cubic splines
C00087 00008	procedure filldraw(real pensize)
C00091 00009	The \\{drawit} procedure begins here
C00098 00010	the (temporary) routine for binary input
C00100 00011	Now comes very system-dependent code for displaying on the user's screen
C00106 ENDMK
C⊗;
entry; begin comment The semantics module of METAFONT.

(It is wise to be familiar with the memory management routines of MFSYS
before going very deeply into the following code.)

This module contains routines that draw characters on a bit raster. It is
pretty much independent of the other modules of METAFONT, although you should
read the documentation of \\{scanpath} in MFNTRP if you want to understand what
information is passed to the \\{drawit} routine. 

The algorithms have been designed to take time (and space) of order $n↑2$
if the resolution is $n$ pixels per inch. The space and time requirements have
also been reduced by using full-word boolean operations. Thus, it is hoped
that METAFONT will be sufficiently quick on the draw;

comment Certain bits of the "control" variable govern the on-line output:
	'100	trace numeric data of curves to be drawn
	'200	trace changes made to avoid sharp corners or high speeds
	'400	trace the points plotted
    '2000000	display the raster after every "draw" or "ddraw"

Furthermore the "circlemode" bit changes the kinds of pen generated:
'10000000000	ddrawn circles replace ellipses
;

require "MFHDR.SAI" source_file;
define drawtrace=⊂(control land '100)⊃, modtrace=⊂(control land '200)⊃,
	plottrace=⊂(control land '400)⊃;
define drawdisplay=⊂(control land '2000000)⊃;
define displaymodes='3000000;
define circlemode=⊂control land '10000000000⊃;
internaldef all_ones=-1 # full word of one bits;

internaldef xxtr=⊂realparam[1]⊃, xytr=⊂realparam[2]⊃, xtr=⊂realparam[3]⊃,
	yxtr=⊂realparam[4]⊃, yytr=⊂realparam[5]⊃, ytr=⊂realparam[6]⊃;
comment These are the coefficients of transformation: we plot
	(xxtr.x + xytr.y + xtr, yxtr.x + yytr.y + ytr)
instead of $(x,y)$;
comment Procedures for generating and storing pens;

internaldef hw=bitsperwd div 2 # number of bits in a halfword;
define lefthalf(x)=⊂(x ash -hw)⊃ # left halfword;
define righthalf(x)=⊂(x land ((1 lsh hw)-1))⊃ # right halfword;

comment The major activity of METAFONT's raster module is to take digitized "pens"
(or erasers) and to make digital images by placing the pen in some specified
point and setting the corresponding bits of the picture to 1 (or resetting them
to 0).

The digitized pen images are stored in a rather peculiar fashion, in order
to facilitate "dynamic" plotting, since the main application is to plot a pen that
is making a sequence of king moves (one step horizontally, vertically, or
diagonally). Eight bit patterns are stored for each pen, corresponding to the
eight possible moves. If $P(x,y)$ denotes the precidate that the pen includes
point $(x,y)$, the eight edge patterns are

	N	P(x,y) ∧ ¬P(x,y+1)
	NE	P(x,y) ∧ ¬P(x+1,y+1)
	E	P(x,y) ∧ ¬P(x+1,y)
	SE	P(x,y) ∧ ¬P(x+1,y-1)
	S	P(x,y) ∧ ¬P(x,y-1)
	SW	P(x,y) ∧ ¬P(x-1,y-1)
	W	P(x,y) ∧ ¬P(x-1,y)
	NW	P(x,y) ∧ ¬P(x-1,y+1)

Thus, for example, the edge patterns corresponding to a pen that is a $3\times3$
square array would be

	 N     NE    E     SE    S     SW    W     NW

	111   111   001   001   000   100   100   111
	000   001   001   001   000   100   100   100
	000   001   001   111   111   111   100   100

This information is redundant, essentially four times as much as needed, since any
of the pairs N-S, NE-SW, E-W, NW-SE gives enough information to deduce the entire
pen shape. But all eight patterns are stored, to gain speed.

For simplicity it is assumed that the pen shapes are "horizontally convex", i.e.,
that all the 1's in each row are consecutive. Pen images are usually generated in
the arrays \\{penl} and \\{penr}, specifying the leftmost and rightmost
bit positions of the pen in each row. More precisely, the integer variables \\{ymin}
and \\{ymax} and the integer arrays \\{penl} and \\{penr} are set so that, for
$\\{ymin}≤y≤\\{ymax}$, the pen image includes points $(x,y)$ if and only if
$\\{penl}[y]≤x≤\\{penr}[y]$. For example, a $3\times3$ pen would be represented
by $\\{ymin}=-1$, $\\{ymax}=+1$, $\\{penl}[y]=-1$ and $\\{penr}[y]=+1$ for
$-1≤y≤+1$.

After a pen has been generated, its eight representations are stored into a big
array called \\{pmem}. This array is allocated sequentially, and it is used to
save all pen images except the special ones defined by \&{spen} or \&{epen}.
Within this array the pens of a given type are linked together in a
sorted, doubly linked list. For example, all pens of type \\{hpen} are
accessible from location \\{hpenhead}, which contains a forward pointer to
the smallest \\{hpen}, which in turn contains both a backward pointer to
\\{hpenhead} and a forward pointer to the second-smallest \\{hpen}, and so on.
The largest \\{hpen} has a forward pointer to $\\{hpenhead}+2$, where the list ends.
In \\{pmem} these forward and backward pointers are followed by the pen size and
a specification of the individual edges.

More precisely, a digitized pen is represented in \\{pmem} as follows:

		lefthalf		righthalf

0th word:	backward pointer	forward pointer
1st word: the pen size, a positive integer less than \\{infty}
2nd word:	   N			   NE
3rd word:	   E			   SE
4th word:	   S			   SW
5th word:	   W			   NW

Here N, NE, etc., are pointers to the first \\{pmem} entries for the corresponding
edge patterns. Each edge pattern (except for the E and W edges) is stored as a
packed bit array in the following form: First comes a word containing
$\\{lefthalf}=\\{xdelta}$, $\\{righthalf}=\\{xcount}$. The bits are stored in
columns of words, with \\{bitsperwd} bits in each word. For example, if each word
contains 36 bits, column zero contains bits for $-17≤x≤+18$, and column $-1$
contains bits for $-53≤x≤-18$, etc. If \\{xlc} and \\{xrc} denote the leftmost
and rightmost column numbers of the edge pattern, then $\\{xdelta}=\\{xlc}$ and
$\\{xcount}=\\{xrc}-\\{xlc}$.\xskip(Note that \\{xdelta} can be negative but
\\{xcount} is always $≥0$. This property will hold in general when two
quantities are being packed into halfwords in METAFONT's data structures.)\xskip
The first word of the edge representation is followed by $\\{xcount}+1$ groups
of words specifying the individual columns from left to right as follows: The
first word of a column group contains $\\{lefthalf}=\\{ydelta}$ and $\\{righthalf}
=\\{ycount}$, where \\{ydelta} is the $y$ coordinate of the bottom word in the
column and $\\{ydelta}+\\{ycount}$ is the $y$ coordinate of the top word. The
following $\\{ycount}+1$ words contain the actual full-word bit patterns of the
pen edge from bottom to top of that column. (If the column is blank, \\{ycount}
will be zero and the following word contains all zeroes.)

For example, suppose we are representing columns with only two bits per word.
Then the NW edge of the $3\times3$ pen discussed above would be presented as a total
of seven words:

	xdelta=-1	xcount=1
	ydelta=-1	ycount=2
		   01
		   01
		   01
	ydelta=1	ycount=0
		   11

The E and W edges are stored in another form intended to facilitate horizontal
painting. The E-W form is simply a word containing $(\\{ydelta},\\{ycount}$,
followed by $\\{ycount}+1$ words containing the (unique) $x$ coordinate of the
pen edge in row $y$, for $y=\\{ydelta}$, $\\{ydelta}+1$, $\ldotss$, $\\{ydelta}
+\\{ycount}$. For example, the E pattern of the $3\times3$ pen would have
$\\{ydelta}=-1$ and $\\{ycount}=2$, with the sequence of $x$ coordinates $+1$,
$+1$, $+1$.
;

internaldef xpenmin=-161,xpenmax=126,ypenmax=99,ypenmin=-ypenmax # legal pen range;
comment $\\{xpenmin}-1$ and \\{xpenmax} should be congruent to \\{hw},
	modulo \\{bitsperwd};
define column(x)=⊂((x-xpenmin) div bitsperwd)⊃;
define bitloc(x)=⊂((x+(1000*bitsperwd+hw-1))mod bitsperwd)⊃ #
	number of bits to the left of bit \\x in a packed word;
saf integer array penl,penr[ypenmin:ypenmax] # pens generated here;
integer ymin, ymax # lower and upper boundaries of a generated pen;
integer xmin, xmax # left and right boundaries of a generated oval pen;

internaldef pmemsize=30000 # number of words of pen storage;
saf integer array pmem[0:pmemsize+8] # storage of pens;
internal integer pmemptr # pointer to first available place in \\{pmem};
integer pmemplace # value of \\{pmemptr} on most recent entry to \\{storepen};
saf integer array penlist[0:column(xpenmax)] # array used by \\{storepen};
integer biggestpen # largest amount of pmem per pen, so far;

define N_dir=⊂lefthalf(pmem[curploc+2])⊃, NE_dir=⊂righthalf(pmem[curploc+2])⊃,
	E_dir=⊂lefthalf(pmem[curploc+3])⊃, SE_dir=⊂righthalf(pmem[curploc+3])⊃,
	S_dir=⊂lefthalf(pmem[curploc+4])⊃, SW_dir=⊂righthalf(pmem[curploc+4])⊃,
	W_dir=⊂lefthalf(pmem[curploc+5])⊃, NW_dir=⊂righthalf(pmem[curploc+5])⊃;
	
procedure storepen # takes generated pen and moves it into \\{pmem};
begin comment Assuming that a nonempty pen is described by \\{ymin}, \\{ymax},
\\{penl}, and \\{penr} as explained above, this procedure converts it to the format
required in \\{pmem}. During the calculation a number of linked lists are used:
\\{penloc}[\\{xc}] points to a list of two-word entries for column \\{xc}, where
node $p$ contains three fields:
	vmemint(p) is a nonzero bit pattern,
	info(p) is the $y$ coordinate minus \\ypenmin,
	link(p) points to the next node in the list.
These lists are in increasing order by $y$ coordinate;
integer i,j # increments to $x,y$ for the current edge;

if pmemptr+2*biggestpen>pmemsize then clearpens(false);
pmemplace←pmemptr;
pmemptr←pmemptr+6 # move pmemptr to the location of the first edge;
for i←-1 thru 1 do for j←-1 thru 1 do
	begin case 3*i+j+4 of begin comment (This +4 is here only because
		SAIL doesn't allow negative cases.);
	[0] pmem[pmemplace+4]←pmemptr # SW;
	[1] pmem[pmemplace+5]←pmemptr lsh hw # W;
	[2] pmem[pmemplace+5]←pmem[pmemplace+5]+pmemptr # NW;
	[3] pmem[pmemplace+4]←pmem[pmemplace+4]+(pmemptr lsh hw) # S;
	[5] pmem[pmemplace+2]←pmemptr lsh hw # N;
	[6] pmem[pmemplace+3]←pmemptr # SE;
	[7] pmem[pmemplace+3]←pmem[pmemplace+3]+(pmemptr lsh hw) # E;
	[8] pmem[pmemplace+2]←pmem[pmemplace+2]+pmemptr # NE;
	else continue
	  end;
	if j then
		begin comment Storing the packed bit image (not E-W) form;
		integer xminc,xmaxc # smallest and largest columns;
		integer xc # runs thru the columns;
		integer y # runs thru the rows;
		xminc←infty;
		for y←ymax step -1 until ymin do
			begin integer l,r,ll,rr; label one_interval,two_intervals;
			l←penl[y]; r←penr[y];
			if (y=ymin and j<0) or (y=ymax and j>0) then
				go to one_interval;
			ll←penl[y+j]-i; rr←penr[y+j]-i;
			comment We want to represent all bits $x$ such that
				$l≤x≤r$ and not $\\{ll}≤x≤\\{rr}$;
			if ll>r or rr<l then go to one_interval;
			if ll≤l then
				begin if rr≥r then continue;
				l←rr+1; go to one_interval;
				end
			else if rr≥r then
				begin r←ll-1; go to one_interval;
				end;
			two_intervals: comment $l≤x<\\{ll}$ and $\\{rr}<x≤r$;
			define interval(l,r)=⊂begin integer xlc,xrc,p,bl,yy;
			xlc←column(l); xrc←column((r)); yy←y-ypenmin;
			if xlc<xminc then
				begin if xminc=infty then
					begin xminc←xmaxc←xlc; penlist[xlc]←0;
					end
				else while xlc<xminc do
					begin xminc←xminc-1; penlist[xminc]←0;
					end;
				end;
			while xrc>xmaxc do
				begin xmaxc←xmaxc+1; penlist[xmaxc]←0;
				end;
			bl←bitloc(l);
			while xlc<xrc do
				begin getvavail(p); vmemint(p)←all_ones lsh(-bl);
				mem[p]←penlist[xlc]+(yy lsh infod); penlist[xlc]←p;
				xlc←xlc+1; bl←0;
				end;
			getvavail(p); vmemint(p)←(all_ones lsh((bitsperwd-1)
				-bitloc((r))+bl)) lsh(-bl);
			mem[p]←penlist[xlc]+(yy lsh infod); penlist[xlc]←p;
			end⊃ # end of the definition of interval(l,r);
			interval(l,ll-1); l←rr+1;
			one_interval: comment $l≤x≤r$; interval(l,r);
			end;
		pmem[pmemptr]←((xminc-column(0)) lsh hw) + xmaxc-xminc;
		pmemptr←pmemptr+1;
		for xc←xminc thru xmaxc do if penlist[xc] then
			begin integer p,pmemp,y,q;
			p←penlist[xc]; y←info(p);
			pmem[pmemptr]←((y+ypenmin) lsh hw)-y;
			pmemp←pmemptr; pmemptr←pmemptr+1;
			loop	begin comment Now $y=\info(p)$;
				label advance_p # go here when finished with $p$;
				pmem[pmemptr]←vmemint(p);
				pmemptr←pmemptr+1;
				if pmemptr≥pmemsize then overflow(pmemsize);
				advance_p: q←p; p←link(p); freeavail(q);
				if p=0 then done;
				if y=info(p) then
					begin comment Two entries for the same $y$;
					pmem[pmemptr-1] ←
						pmem[pmemptr-1] lor vmemint(p);
					go to advance_p;
					end;
				loop	begin y←y+1;
					if y=info(p) then done;
					pmem[pmemptr]←0;
					pmemptr←pmemptr+1;
					if pmemptr≥pmemsize then overflow(pmemsize);
					end;
				end;
			pmem[pmemp]←pmem[pmemp]+y;
			end
		else	begin comment Empty list;
			pmem[pmemptr]←pmem[pmemptr+1]←0; pmemptr←pmemptr+2;
			if pmemptr≥pmemsize then overflow(pmemsize);
			end;
		end
	else	begin comment The simple E-W form;
		integer y # runs through the rows;
		pmem[pmemptr]←(ymin lsh hw)+ymax-ymin;
		pmemptr←pmemptr+1;
		if pmemptr+ymax-ymin≥pmemsize then overflow(pmemsize);
		for y←ymin thru ymax do
			begin pmem[pmemptr]←(if i<0 then penl[y] else penr[y]);
			pmemptr←pmemptr+1;
			end;
		end;
	end;
if pmemptr-pmemplace>biggestpen then biggestpen←pmemptr-pmemplace;
DEBUGONLY print(" [",pmemptr,"]");
end;

procedure makeovalpen(real a,b,c,x0,y0;boolean nostore)
	# generate a pen with elliptical shape;
begin comment This procedure creates a pen containing all integer points $(x,y)$
such that $a(x-x0)↑2+b(x-x0)(y-y0)+c(y-y0)↑2≤1$, assuming that $a>0$, $c>0$,
and $b↑2-4ac<0$. If nostore is false, the pen is stored in the \\{pmem} array
beginning at the location stored in \\{pmemplace};
real delta # the discriminant;
real s # maximum of $|y-y0|$;
real t # maximum of $|x-x0|$;
integer y # runs through possible \\y values;
label bad # go here when the assumptions are violated;
label ready # go here when a pen has been generated ready to be stored;
if (delta←4*a*c-b↑2)≤0.0 then go to bad;
if a<0 then go to bad; if c<0 then go to bad;
s←sqrt(4*a/delta); t←sqrt(4*c/delta);
ymax←floor(s+y0); ymin←-floor(s-y0);
if ymax>ypenmax then overflow(ypenmax);
if ymin<ypenmin then overflow(ypenmin);
xmin←infty; xmax←-infty;
for y←ymin thru ymax do
	begin real d,yy; integer xl,xr;
	yy←y-y0;
	d←4*a-delta*yy↑2; if d<0 then d←0 else d←sqrt(d);
	xr←penr[y]←floor((-b*yy+d)/(2*a)+x0);
	xl←penl[y]←-floor((b*yy+d)/(2*a)-x0);
	if xl<xmin and xl≤xr then xmin←xl;
	if xr>xmax and xl≤xr then xmax←xr;
	end;
if xmin=infty then go to bad # empty pen;
while penl[ymin]>penr[ymin] do ymin←ymin+1;
while penl[ymax]>penr[ymax] do ymax←ymax-1;
if xmin<xpenmin then overflow(xpenmin);
if xmax>xpenmax then overflow(xpenmax);
go to ready;
bad: error("Empty pen specification");
penl[0]←penr[0]←ymin←ymax←xmin←xmax←0 # substitute single bit at (0,0);
ready: if not nostore then storepen;
end;

procedure makepen(integer xwidth,ywidth) # generate a standard oval pen;
begin comment This procedure makes an elliptical pen of dimensions
$\\{xwidth}\times\\{ywidth}$ pixels, assuming that both parameters
are positive integers;
integer x,y;
real xcorr,ycorr # pen offsets;
real xw,yw,xc,yc # parameters sent to makeovalpen;
real fudge # A factor that guarantees correct width and height in terms of pixels;
x←(xwidth+1) land 1; if x then xcorr←.5 else xcorr←0.0;
y←(ywidth+1) land 1; if y then ycorr←.5 else ycorr←0.0;
xw←xwidth; yw←ywidth; xc←xcorr; yc←ycorr;
if circlemode then
	if xwidth≥ywidth then
		begin xw←ywidth; x←y; xc←ycorr;
		end
	else	begin yw←xwidth; y←x; yc←xcorr;
		end;
fudge←x/xw; if y/yw>fudge then fudge←y/yw;
fudge←1.0/(1.0+fudge↑2);
makeovalpen(fudge*(2.0/xw)↑2, 0, fudge*(2.0/yw)↑2,xc,yc,circlemode);
if circlemode then
	begin if xwidth>ywidth then
		begin comment stretch the circle horizontally;
		integer sl,sr,k;
		sl←(xwidth-ywidth)/2+ycorr-xcorr+.2;
		sr←xwidth-ywidth-sl;
		xmin←xmin-sl; xmax←xmax+sr;
		if xmin<xpenmin then overflow(xpenmin);
		if xmax>xpenmax then overflow(xpenmax);
		for k←ymin thru ymax do
			begin penl[k]←penl[k]-sl; penr[k]←penr[k]+sr;
			end;
		end
	else if xwidth<ywidth then
		begin comment stretch the circle vertically;
		integer sd,su,k;
		sd←(ywidth-xwidth)/2+xcorr-ycorr+.2;
		su←ywidth-xwidth-sd;
		ymin←ymin-sd; ymax←ymax+su;
		if ymin<ypenmin then overflow(ypenmin);
		if ymax>ypenmax then overflow(ypenmax);
		for k←ymin step 1 until -1-sd do
			begin penl[k]←penl[k+sd]; penr[k]←penr[k+sd];
			end;
		for k←ymax step -1 until 1+su do
			begin penl[k]←penl[k-su]; penr[k]←penr[k-su];
			end;
		for k←-sd thru su do
			begin penl[k]←penl[0]; penr[k]←penr[0];
			end;
		end;
	storepen;
	end;
end;

procedure makerpen(integer xwidth,ywidth) # generate a standard rectangular pen;
begin comment This procedure makes a rectangular pen of dimensions
$\\{xwidth}\times\\{ywidth}$ pixels, but shifted to the left or right of the
origin depending on whether \\{curpen} is \\{lpen} or \\{rpen},
assuming that both parameters are positive integers;
integer xl,xr,y;
if xwidth>-xpenmin then
	begin error("Rectangle too wide"); xwidth←-xpenmin;
	end;
xl←-xwidth; xr←xwidth;
ymin←1+((-ywidth) ash -1); ymax←ywidth ash -1;
if ymin<ypenmin then overflow(ypenmin);
if curpen=lpen then xr←-1 else xl←1;
for y←ymin thru ymax do
	begin penl[y]←xl; penr[y]←xr;
	end;
storepen;
end;

internal saf real array spenspec[1:7] # specifications for a special pen;
real sxcorr,sycorr,sxmin,sxmax,symin,symax # pen offsets for special pen;

internal procedure makespen # create a new special pen;
begin makeovalpen(spenspec[1],spenspec[2],spenspec[3],spenspec[4],spenspec[5],false);
sxcorr←spenspec[6]; sycorr←spenspec[7];
sxmin←xmin; sxmax←xmax; symin←ymin; symax←ymax;
pmemptr←pmemplace # this position in \\{pmem} will be reusable;
end;
	
internaldef epensize=ypenmax-ypenmin+1 # maximum length of \&{epen} specs;
internal saf integer array epenlspec,epenrspec[0:epensize] # explicit pen specs;
internal integer epen0,epenptr # zero point and end of explicit pen specs;
real exmax,exmin,eymax,eymin # pen offsets for explicit pen;

internal procedure makeepen # create a new explicit pen;
begin integer k,m,y;
k←epen0-1; m←epenptr-epen0;
if epenyfactor≤0 then
	begin error("epenyfactor must be positive (1.0 assumed)"); epenyfactor←1.0;
	end;
if epenxfactor≤0 then
	begin error("epenxfactor must be positive (1.0 assumed)"); epenxfactor←1.0;
	end;
ymax←epenyfactor*k;
ymin←epenyfactor*m; ymin←-ymin # truncate towards zero;
if ymax>ypenmax then overflow(ypenmax);
if ymin<ypenmin then overflow(ypenmin);
xmin←infty; xmax←-infty;
for y←ymin thru ymax do
	begin real yy,alpha; integer t; yy←y/epenyfactor;
	t←yy; alpha←yy-t; t←epen0-t;
	penl[y]←epenxfactor*(epenlspec[t]+alpha*(epenlspec[t-1]-epenlspec[t]))+.5;
	penr[y]←epenxfactor*(epenrspec[t]+alpha*(epenrspec[t-1]-epenrspec[t]))+.5;
	if penl[y]<xmin then xmin←penl[y];
	if penr[y]>xmax then xmax←penr[y];
	end;
if xmax>xpenmax then overflow(xpenmax);
if xmin<xpenmin then overflow(xpenmin);
storepen;
exmin←xmin; exmax←xmax; eymin←ymin; eymax←ymax;
pmemptr←pmemplace # this position in \\{pmem} will be reusable;
end;

internal integer curploc # location of the current pen in \\{pmem};
internal boolean eraser # the current "pen" really is an eraser;
comment When the current pen type (\\{curpen}) changes, \\{curploc} should
be set to zero, and \\{eraser} should be set to \false. The routines assume that,
if $\\{curploc}≠0$, it points to a pen of type \\{curpen} but possibly of the
wrong width, and that \\{eraser} is properly set;

define cpenhead=0,hpenhead=4,vpenhead=8,lpenhead=12,rpenhead=16;
define pmemstart=20 # first location in \\{pmem} that isn't preallocated;

internal procedure resetspen # forgets the previous spen specification, if any;
spenspec[1]←spenspec[2]←spenspec[3]←spenspec[4]←spenspec[5]←spenspec[6]←
	spenspec[7]←0;

internal procedure resetepen # forgets the previous epen specification, if any;
begin epenlspec[1]←epenrspec[1]←0; epen0←epenptr←1;
end;

internal procedure resetpens # the current pen is unknown;
begin curpen←badpen; cursize←0.0; curploc←0; eraser←false;
end;

internal procedure clearpens(boolean all) # initializes the pen memory;
begin comment This procedure wipes out all the information currently in \\{pmem},
and does all other initialization required for the pen/eraser routines;
integer h # runs through cpenhead, hpenhead, etc.;
if all then
	begin resetpens; resetspen; resetepen; biggestpen←0;
	end
else curploc←0;
for h←cpenhead step 4 until pmemstart-4 do
	begin comment Set doubly linked lists empty;
	pmem[h]←h+2; pmem[h+1]←0;
	pmem[h+2]←h lsh hw; pmem[h+3]←infty;
	end;
pmemptr←pmemstart;
end;

procedure setuppen(integer w) # ensures that the current pen is in \\{pmem};
begin comment This procedure makes sure that \\{curploc} is pointing to the
pen node of type \\{curpen} and size \\w;
if curploc=0 then
	begin case curpen of begin
	[cpen] curploc←cpenhead;
	[hpen] curploc←hpenhead;
	[vpen] curploc←vpenhead;
	[lpen] curploc←lpenhead;
	[rpen] curploc←rpenhead;
	[spen] begin makespen; curploc←pmemplace end;
	[epen] begin makeepen; curploc←pmemplace end;
	else begin error("No pen defined"); curpen←cpen # $\\{cpenhead}=0$; end
	  end;
	if ((w≤0) or (w≥infty)) then
		begin error("Illegal pen size ("&cvs(w)&")"); w←1;
		end;
	end;
if curpen≥spen then
	begin comment The size of an spen or epen is ignored;
	pmem[curploc+1]←w; return;
	end;
if pmem[curploc+1]=w then return;
if plottrace then print("|",w,"|");
if pmem[curploc+1]<w then 
	do curploc←righthalf(pmem[curploc]) until pmem[curploc+1]≥w
else	do curploc←lefthalf(pmem[curploc]) until pmem[curploc+1]≤w;
if pmem[curploc+1]=w then return;

comment It is necessary to generate a new pen image;
case curpen of begin
[cpen] makepen(w,w);
[hpen] makepen(w,hpenht);
[vpen] makepen(vpenwd,w);
[lpen] makerpen(w,lpenht);
[rpen] makerpen(w,rpenht);
else confusion
  end;
if pmem[curploc+1]>w then
	begin comment Insert to left in doubly linked list;
	integer p; p←lefthalf(pmem[curploc]);
	pmem[p]←(pmem[p] land (all_ones lsh hw))+pmemplace;
	pmem[curploc]←(pmem[curploc] land (all_ones lsh-hw))+(pmemplace lsh hw);
	pmem[pmemplace]←(p lsh hw) + curploc;
	end
else	begin comment Insert to right in doubly linked list;
	integer p; p←righthalf(pmem[curploc]);
	pmem[p]←(pmem[p] land (all_ones lsh-hw))+(pmemplace lsh hw);
	pmem[curploc]←(pmem[curploc] land (all_ones lsh hw))+pmemplace;
	pmem[pmemplace]←(curploc lsh hw) + p;
	end;
curploc←pmemplace; pmem[pmemplace+1]←w;
end;

internal real procedure penadj(real width; integer dir) # boundary of pen position;
begin comment For example, "lft9 z" equals $z$ plus $\\{penadj}(w↓9,\.{lft})$;
integer w;
w←width+.5; if w≤0 then w←1;
if curpen≥spen then setuppen(1) # ensure that an spen or epen has been generated;
case 4*curpen+dir of begin
[4*cpen+lft][4*cpen+bot][4*hpen+lft][4*vpen+bot] return(.5*(1-w));
[4*cpen+rt][4*cpen+top][4*hpen+rt][4*vpen+top] return(.5*(w-1));
[4*hpen+bot] return(.5*(1-hpenht));
[4*lpen+bot] return(.5*(1-lpenht));
[4*rpen+bot] return(.5*(1-rpenht));
[4*vpen+lft] return(.5*(1-vpenwd));
[4*hpen+top] return(.5*(hpenht-1));
[4*lpen+top] return(.5*(lpenht-1));
[4*rpen+top] return(.5*(rpenht-1));
[4*vpen+rt] return(.5*(vpenwd-1));
[4*lpen+lft] return(-w);
[4*lpen+rt] return(-1.0);
[4*rpen+lft] return(1.0);
[4*rpen+rt] return(w);
[4*spen+lft] return(sxmin-sxcorr);
[4*spen+rt] return(sxmax-sxcorr);
[4*spen+bot] return(symin-sycorr);
[4*spen+top] return(symax-sycorr);
[4*epen+lft] return(exmin-excorr);
[4*epen+rt] return(exmax-excorr);
[4*epen+bot] return(eymin-eycorr);
[4*epen+top] return(eymax-eycorr);
else comment do nothing;
  end;
error("Undefined pen"); return(0.0);
end;

string procedure strpen # symbolic name of current pen;
return(case curpen of("cpen","hpen","vpen","lpen","rpen","spen","epen","no pen"));
comment Procedures for plotting points on the raster;

comment METAFONT draws characters in a big array called \\{rast}, containing
bits packed into words. When a pen is positioned at point $(x,y)$, its
pixel coordinates are shifted by addition of $x$ and $y$, to obtain the
points of \\{rast} that are set or reset (depending on whether \\{eraser} is
\\{false} or \\{true}. Only points $(x,y)$ in the range $\\{xrastmin}≤x
≤\\{xrastmax}$ and $\\{yrastmin}≤y≤\\{yrastmax}$ may be plotted.
Here \\{xrastmin} and \\{xrastmax} are multiples of \\{bitsperwd}.
Bit $(x,y)$ in the raster appears in word $\\{rast}[\\{rloc}(x,y)],
with $\\{bitloc}(x)$ bits to its left in that word.
The allocation is such that bits $(x,\\{yrastmax}+1)$ thru
$(x,\\{yrastmax}+\\{ypenmax})$ occupy the same position as bits
$(x+\\{bitsperwd},\\{yrastmin}+\\{ypenmin})$ thru
$(x+\\{bitsperwd},\\{yrastmin}-1)$;

comment No procedures except those on this page are supposed to change \\{rast}.
Note that the rast array is arranged very counterintuitively.  Words stretch
out horizontally to form rows, as one might expect, but consecutive words form
columns.  Thus to move to the next word in a row, add rspan to your current
word address.  

The routines maintain four variables \\{xleft}, \\{xright}, \\{ylow}, and
\\{yhigh} such that $\\{rast}[k*rspan+l]$ is zero unless
	xleft ≤ k ≤ xright   and   ylow ≤ l ≤ yhigh;

internaldef xrastmin=-72,xrastmax=360,yrastmin=-75,yrastmax=300 # raster bounds;
internaldef rspan=yrastmax+ypenmax+1-yrastmin # words per raster column;
internaldef rcol(x)=⊂((x-(xrastmin+xpenmin)) div bitsperwd)⊃ # column for bit $x$;
internaldef rloc(x,y)=⊂rcol(x)*rspan+y⊃ # allocation function for \\{rast};
internaldef rast0=rloc(xrastmin+xpenmin,yrastmin+ypenmin),
	rast1=rloc(xrastmax+xpenmax,yrastmax+ypenmax)+rspan;
require nextline&"The raster contains "&cvs(rast1+1-rast0)&" words"&nextline message;
comment SAIL won't let the above exceed 2↑17=131072;
comment At PARC, due to a TENEX Sail restriction, the entire dirty part
	of your core image can't exceed 128K, so the raster must be
	substantially smaller.  In addition, the Loader doesn't work if
	the entire core image being loaded exceeds 128K.  Therefore, at
	PARC, the raster array must be allocated dynamically.  On the
	other hand, we also want to declare it in the outermost block.
	The way around this seeming dilemma is to make the raster array a
	field of a record that is dynamically allocated.  This makes
	the code run a little slower, but keeps the Loader happy.
	To understand the following hackery, look at page 65 of the
	Sail manual;
IFWAITS
internal saf integer array rast[rast0:rast1] # the big raster workspace;
ENDWAITS
IFTENEX
internal record_class rastrec(saf integer array r);
internal record_pointer(rastrec)rastptr;
internaldef rast = ⊂rastrec:r[rastptr]⊃;
ENDTENEX
IFTOPS20
IFXMEM comment The raster will live in sections 2 and up.  Look at the last
	 page of this file for all the macros used to get to it;
ELSEC
internal saf integer array rast[rast0:rast1] # the big raster workspace;
ENDC
ENDTOPS20

internal integer xleft,xright,ylow,yhigh # active part of wkgO&+Il4Ph+≠?↔;πK⊃πβK?∂.#WK∃ε≠3↔π⊗#⊃↓
ε≠3↔π∩βS#∃∧#πSπ&KO
β↔+≠≠↔∪X4+'w#↔K;∞aβCK}≠↔∪W⊗)β∂3.KKπ∨!↓
β≡+SMβ⊗OS↔∩βS=βV+K=lhS↔∨Nqβ';&+∨↔IεY31lhS≠?IεZ␈c3.3QβSG∪Uβc⊗K∨#Qε#<4(HK↔∨Nqβ';&+∨↔Iεc?
mεc?∞␈ZSKOC∞q/g3␈9l4(HKKπO"∨↔S~↔cC∩C3?
c↓%l4TJ~b6,h$'c⊗cS'QG∪πOQ∞c?
#f{
-EJcKπO"3?
Fc?
%gK#'∨Bkg3?:Il4*,bN⊗HKπKK⊗cQ#K∂≠Ro3}→-Fug∪πOR↑c?∞ugK#'∨Bkg3?:Il4*,r∩4PH'∂?nk↔;Qε3?Iβe{g3?:βS#K*βg#'>Aβ∪=h($$O∪πOQ∞;↔SM∞+cCIFY+KOε9/1c↓%↓
π∪πOR↑Y+KOε9/2myAl4PH'↔;#X4+cf+≠R␈Nc?↑␈Ns≠SeZβcK'>CR␈gFK∨"⎇nK;≠SKX4+∂f+πK∪#X4+↔v!l4(hSCK?≡+∪WK*β[Kπ∨#C3?"C';S.;↔IβαcaA3K↓39%α→βC↔r↓∪A⊃εCC3N+⊃βSz↓#aAgIA%9rq#aAgIA/9KX4+.;'9β≡{77↔w!αS#O→βCK}≠↔∪W⊗)βWO/→βS#*βC↔9ε{Iβ↔⊗O↔Iπ≠C↔∂N3'↔⊃ε∪eβSF)β;?rj∃6\hS↔∪∨*βCπS&+K9β⊗+∨';vK;≥β∂!β3?≡S'?r↓∪A⊃εK9αrg[C7↔oqβS=π+C∪π&)βS#*βKπO&+H4+∂!↓∪9[	⊃βC}K;SMεK9β¬π3↔KSN≠π1βfK;∃β>{';≥π+Aβ≠⊗{5βC}K;Q↓"CaA3K↓%⊃9∧KQβπ∨≠W7↔_h+S#∂!βπ3bβ?→β&C↔O∃πβ?';'→βπK*β'9β⊗;∨∃αCS#∃ε≠π33Ns≥βK␈+S';*βO#?.c⊃β∂F+∂-β&C'M%Xh+';&+∨↔IπC	3c:cc∂?.sQ3cf→3cK~cg13NAl4+N1βC3␈#SKπ≡)βS#.qβCKNsQ!	B⊃3aAb⊃1	3K↓1	i∩ceA/ra	%	KX4+c∃y#aAoCKπO&k'9'n{⊃βO#OC↔↔;⊃mβFc∞␈K≡{1#aαkc	%↑c↔≠SF3→#εk↔6o¬i%l4WC∂?Ww"␈K'>CS#πf1#C7.joBuJ↓
βK∞s∨∃β}1β∂?g+7;Mπ≠Cπ;v+⊃βJβS#∃πβ↔9lhScK∞␈C3
/F≠?W;#X4+'2βc3
gC3↔≠"βS#↔rβc3↔7"␈c3≠X4+'2βcKoCK'∨G!βS#.qβcKN;#R␈G∪
-EXh+c↑␈C3
+↔≠Cπ9α→β3?≡S'?rβ?→βε+9β'rαrso⊗OSyXh+B␈αYEl4Vc??@N∪↔∨'rβ';S.;↔IβGI↓
β∨+KK↔w!β←?⊗!β'9¬bsoK∂≠SylhP'';&+∨↔IεY↓
αecog∂␈+;SyεK9β∂/∪K↔;"β∂?3.k9l4PK';S.;↔IβJc)↓
εc??AεK;∪'≡+M1β↔+;;'v9βS#↔)βC↔rβ∂?3.k;MlhP'';&+∨↔IπQ↓
β∨+KK↔w!βC?↔#'?9ε{→βC.qβ'7∞;∃l4PKg2␈K↓/3↔7##π32CC7↔m[Bu%Zβ.␈KN;#S#∞c→#Cn+6oBjIl4(OK"␈gb[-/9Xh(''2βg1sNc?]β&C↔9βNc?↑␈Nal4(NK→βgCsg#'>AβS#.qβg#N;"␈gCX4('GJ␈c]←K1mβUyAl4PK∂?7n+;Qαv{]β∂}k↔Mβ&C∃ββεK;;↔∩β3??α9≥l4PK∪↔≠Ns∃βππβ3gC.qtC'2β↔Kπ≡+IβSF+9β.;'84PH'Kπ∨!π∨↔'→πKπ∨!π3πv!π3;␈!π↔cπ⊃#cegQβ3OB↓7c	Hh($%α→βKπ∨"ocfm{KπO%[cfuεcπ;⊃εc;?QαCiβ3≡A↓7c∩Il4(HKKπO"∨↔S~KπO"3π;"3;?"↔cC∩Cce/↔≠Cπ9gQβ3OB↓#''≠C↔K>!7c	JH4($J↓
βK∂≠RocJ[KOC∞rv␈K∂≠RocJ[KOC∞ruβ3∞s⊃β3v{Q↓#Rβ3O!αC'S∨β↔K←"kc	%KX4($OCf␈cJYEβ↔v 4($N+3O∃ε∪↔∨'ph($'⊗OQπ>+SMπ⊗OQπf{Iπ↔GβI#cJciβ3≡A↓7c∩H4($J↓
βK∂≠RocMj␈Kπ∨"ocfjβ3?IαCiβ3≡A↓7c∩Il4(HKKπO"∨↔S~KπO"3?I∞+cCIGCe/K∨βπ93Rβ3O!αC'S∨β↔K←"kc	%Hh($%α→βKπ∨"oce←∪OCπuj␈Kπ∨"oce←∪OCπuiβ3?∩↓#iβg≠!↓#⊗KSOC/∪←⊃7F⊃%%lhP$'cM{ce-
β↔;⊂X4('N1β9sZβS#↔ph($'⊗+∨'9ε3?IβMyAβSG∪Uβ9ε#<4(HH'↔>K9βj␈Qβc?∩βC7↔m[A/%[
umβ∂βC3gε+9l4PH$'↔v!l4(HK≠?IεJ␈9-
βS#K*β-β∪xh($$N∪↔∨'rβj⎇#Rβc?Iπβ7↔6←↓/%-
i%βc␈⊃βC7.joA/Jk:umεCC3Oβ↔9lhP$$'.s⊃l4PH'≠?∩β&␈-nqβS#↔)β-5
β∪<4PH$'.;'9βU{iβc␈⊃βC7.joA/JYFumεCC3Oβ↔9lhP$$'.s⊃l4PH'↔; h('↔g≠∀'.;'9β6{Iβ&{↓βS#↔)β-β&x4($HK↔∨Nqβj␈Rβc?Iπβ7↔6←↓/%-
imβππβ3gC.ql4(HH'↔;#X4($N3?IβMyEβSG∪Uβ9nYβ∪<hP$$'⊗+∨'9εCC3Oβ↔9lhP$$'.s⊃l4PH'≠?∩β&⎇Aπ##KUεY5Eβ&x4($HK↔∨Nqβj␈Rβc?Iπβ7↔6←↓/%-
imβππβ3gC.ql4(HH'↔;#X4($N+;⊃lhP'B␈α[--IXh(''2βc∂?.sQuAπ##↔9ε#?;∃Xh('c≡{W;R␈C∂?Ww!5ElhP'c↑␈C]/K∨βπ9↓~βπ∪[∞s∂∃β&yβ;↔G!β∂?g+79lhP'↔;#X4+↔v!l4(hSCK?≡+∪WK*β#Kπ∨#C3?"C';S.;↔IβαcaA3K↓39%α→βC↔r↓∪A⊃εCC3N+⊃βSz↓#aAgIA%9rq#aA↑q3eAKX4+.;'9β≡{77↔w!αS#O→βCK}≠↔∪W⊗)βWO/→βS#*βC↔9ε{Iβ↔⊗O↔Iε+∪∨∃π≠S?K.!β'9∧)6]β6{K44V∪↔∨'vs';≥εQβ3}≠πS'}q↓∪A"βS=β/β∪πS*βS#∃π∪πOS/⊃βπQα#9-E"βC?'w#Mβ'rβ¬β#␈∪'k?w#π1βfK;∀4V;?';:βS=β&C∃βKN;#Qβ7∪?5βε{';Qα!#aAgIA%⊃rα'Qβ∂≠OW7/→βS#∂!βπ3bβ?→β&C↔O∃πβ?';'→βπK(h+'9π∪π;∨*↓#S#*β∂π3fK;≥β⊗{WS'v)βO#␈+3⊃β≡C↔∂-π##'MKX4+'w#↔∨↔∩βe3g≡{W;QgC13c∩cc3]gCK]3F⊃3g!gC3
3G∪
l4VK→βCf{SSK∞≠∃βSF+9βC⊗K;Q!∩A	3aαa	i	gAA/9b⊃1	3K↓1	%∩Il4+M{eA/f+≠S#∞c→#Cn+6oBjImβg≡{W;R␈∪'∨#&Cπ3→Gβ7↔6←αu%lhSg"␈J[g∂?.sQl4VK→βegK3?]π##↔9πK3?↑␈Il4+N1βg!wK#'∨BβS#↔rβg#'>B␈g!Xh+3?␈'↔>K9βB␈↓-EmπC2␈aα[C7↔m[BumπCJ␈cb[9l4PKc3∞␈∪∂?1GC1%mπCK∞␈⊗≠?1#G⊃%l4PK'→βFc
scf+≠Qβ&C↔9βFc↔≠R␈C3
lhP''→πCK
{G∪'∨#"βS#↔rβcK'>CR␈c⊗→l4(OC3↑␈Fc
+K∨βπ9/KYβcK={cK
W∪OCπr[emβF∩␈'&c?
#Fa%l4PK'→β/∪πO↔∩βS#↔ph($'⊗+∨'9π;#'3*βc3]gCK]β&x4($HK↔∨Nq4(HH'Kπ∨!π∨↔'→πKπ∨!π3πv!π↔cπ⊂4($HH%#cg93π3aC?;↔~β3O!αC'S∨β↔K←"kc	%Hh($$J↓
βK∂≠Rocg:v␈K∂≠Rocg:uβ3∞s⊂4(HH$%#∞c0c?v+Mβ3≡A↓#O#OC↔↔;⊃7c∩I%l4PH$'cg:␈c3:[KOC∞qmβc∃yAl4PH$'↔v!l4(HKKπO"∨↔S~KπO"3π;"↔cC⊂h($$JCc3]bCπ30F{;↔MεcO!↓jC'Sf{
#c∩I7c	[	%'K␈!!7c∩I$4(HI
βK∂≠Rocg:v␈K∂≠Rocg:uβ3∞s⊂4(HH%!#∞c0c?v+Mβ3≡A↓5#⊗KS3?~CcI%oC	-EJKK?QBkc	%KX4($N+;⊂4PK↔3O(K↔∨Nqβ←#Nc∃βcg9scK:β∪<4PH$'.;'9hP$$'⊗OQπ>+SMπ⊗OQπf{Iπ↔GβI#cg93π3aC?;↔~β3O!αkc	$hP$$%~βKπO%[c3↑m{KπO%[c3↑jβ3?IαCπ30F{;↔MεcO!↓oC	%lhP$$'Fc↑␈cg9/KOε9mβF∩⎇AlhP$$'.s⊃l4PH'Kπ∨!π∨↔'→πKπ∨!π3?∩↔cC∩Cc3]`h($$JCπ30F{;↔MεcO!↓BC'S∨β↔K←"iE%7⊗KS3?~CcI%←C	%%εcO!↓oC	$4PH%
β⊗ORoFc↑v␈⊗ORoFc↑uβf{H4(HI↓↓↓αA#π3aC?;↔~β3O!αA#''≠C↔K>!5E%n∪'S3}→#cIJ[c	%Jβ3O!αkc	%Xh($'.s⊃l4PK'→βN≠?W;#iAβSF+9β∪}s∃l4PKg∂?.sR␈g≡{W;Qk	l4(OJ␈e-X4('.s⊃l4V+;⊃lhP4+C⊗{∂↔∪/∪∃β≠.c3Kπ∨#C3?"C';S.;↔IβC↓3eAJ↓
β∂/∪K↔;"βC↔9εQ↓⊃GAA3eαI⊃l4V∪↔∨'rβ∂?7n+;Qα&C'Mβπ∪?∂↔'+K∃β/≠↔Mβ&C∃β∂/∪K↔;"βC↔9ε{)d⊃r eraser to update the
raster at point $(x0,y0)$, including all interior points (not just an edge).
It assumes that point $(x0,y0)$ is in range (the calling routine should check this);
integer p1,p2,xb,xw,xcount,xlc,xrc;
if plottrace then print("(",x0,",",y0,")");
p1←lefthalf(pmem[curploc+4])+2 # S edge of the current pen;
p2←lefthalf(pmem[curploc+2])+2 # N edge of the current pen;
if pmem[p1-2] ≠ pmem[p2-2] then confusion;
xb←(x0-xrastmin)mod bitsperwd; xlc←rcol(x0-xb)+lefthalf(pmem[p1-2]);
xcount←righthalf(pmem[p1-2]) # the number of columns;
xrc←xlc+xcount;
if xlc<xleft then xleft←xlc;
if xrc≥xright then xright←xrc+1;
xw←xlc*rspan # location of pen in \\{rast};
loop	begin integer xy # current word in \\{rast};
	integer q1,q2 # end of current column;
	integer y1,y2 # current row;
	integer z # current portion of pen image;
	define applypen=⊂if eraser then begin
		rast!gets!rast!land!lnot!expr(xy,z lsh -xb)
		# rast[xy]←rast[xy] land lnot (z lsh -xb);
		rast!gets!rast!land!lnot!expr(xy+rspan,z lsh (bitsperwd-xb))
		# rast[xy+rspan]←rast[xy+rspan] land lnot (z lsh (bitsperwd-xb));
		xy←xy+1 end
		else begin
		rast!gets!rast!lor!expr(xy,z lsh -xb)
		# rast[xy]←rast[xy] lor (z lsh -xb);
		rast!gets!rast!lor!expr(xy+rspan,z lsh (bitsperwd-xb))
		# rast[xy+rspan]←rast[xy+rspan] lor (z lsh (bitsperwd-xb));
		xy←xy+1 end⊃;
	y1←lefthalf(pmem[p1-1])+y0; y2←lefthalf(pmem[p2-1])+y0;
	q1←p1+righthalf(pmem[p1-1]); q2←p2+righthalf(pmem[p2-1]);
	xy←xw+y1; z←0;
	if y1<ylow then ylow←y1;
	while y1<y2 do
		begin if p1≤q1 then
			begin z←z xor pmem[p1]; p1←p1+1;
			end;
		applypen; y1←y1+1;
		end;
	while p1≤q1 do
		begin z←z xor pmem[p1]; p1←p1+1;
		applypen; y1←y1+1;
		z←z xor pmem[p2]; p2←p2+1;
		end;
	while p2≤q2 do
		begin applypen; y1←y1+1;
		z←z xor pmem[p2]; p2←p2+1;
		end;
	if z≠0 then confusion;
	if y1-1>yhigh then yhigh←y1-1;
	if xcount=0 then done;
	xcount←xcount-1; xw←xw+rspan;
	p1←p1+1; p2←p2+1;
	end;
end;
comment Beginning of the main procedure \\{drawit};

boolean plotted # we have plotted a point on the current curve;
integer lastx,lasty # if \\{plotted} is \true, these are the coordinates of the
	last point plotted;

saf real array tt[0:9] # boundary values for monoctantic subintervals in \\{spdraw};

define pmin=xrastmin min (-xrastmax) min yrastmin min (-yrastmax),
	pmax=-pmin # legal coordinate range;
saf integer array xl,xr[pmin:pmax] # coordinate bounds computed in \\{spdraw} loop;
saf real array tl,tr[pmin:pmax] # $t$ values computed in \\{spdraw} loop;

saf real array xdel,ydel,del,wdel,tanw[0:maxpoints] # differences used when
	computing spline coefficients;
saf real array w,x,y,xx,yy[0:3,1:maxpoints] # coefficients of cubic splines;

internal procedure drawit(boolean ddrawit) # draws a specified curve;
begin comment This procedure is called on to perform a \.{draw} or \.{ddraw}
operation, after the arrays \\w, \\x, \\y, \\{tanx}, \\{tany}, \\{pointstab},
etc., have been set up as described in \.{MFNTRP}. The code for \\{drawit}
begins several pages below. It is necessary to define several procedures within
\\{drawit} so that errors that abort a drawing will be able to return control
to the main interpretive routine;

label end_of_drawit # go here when drawing is finished or aborted;

procedure drawerror(string s) # give error message and abort the drawing;
begin error(s); go to end_of_drawit;
end;

procedure vplot(integer x0,y0,y1) # current pen applied to $(x0,y0:y1)$;
begin comment This procedure uses the current pen or eraser to update the raster
in a vertical line going up or down from points $(x0,y0)$ to $(x0,y1)$ inclusive;
label plot_one # go here if $(x0,y0)$ must be fully plotted with interior points;
label good_case # go here if the same edge can be used for the whole plot;
label remove_one # go here to remove the first point $(x0,y0)$;
label finished # go here when the plotting is done;
integer p # pointer to proper edge of current pen;
if x0<xrastmin or x0>xrastmax or y0 min y1 < yrastmin or y0 max y1 > yrastmax
	then drawerror("Curve out of range ("&cvs(x0)&","&cvs(y0)&":"&cvs(y1)&")");
if plotted then
	begin integer i,j;
	i←x0-lastx; j←y0-lasty;
	if abs(i)>1 or abs(j)>1 then go to plot_one;
	case 3*i+j+4 of begin
	[0] begin vrastplot(SW_dir,x0,y0,0); go to remove_one end;
	[1] begin hrastplot(W_dir,x0,y0,0); go to remove_one end;
	[2] begin vrastplot(NW_dir,x0,y0,0); go to remove_one end;
	[3] if y1<y0 then go to good_case else begin vrastplot(S_dir,x0,y0,0);
	go to remove_one end;
	[4] go to good_case;
	[5] if y1≥y0 then go to good_case else begin vrastplot(N_dir,x0,y0,0);
	go to remove_one end;
	[6] begin vrastplot(SE_dir,x0,y0,0); go to remove_one end;
	[7] begin hrastplot(E_dir,x0,y0,0); go to remove_one end;
	[8] begin vrastplot(NE_dir,x0,y0,0); go to remove_one end;
	else confusion
	  end;
	end
else plotted←true;
plot_one: fullrastplot(x0,y0);
remove_one: if y0=y1 then go to finished;
if y0<y1 then y0←y0+1
else if y0=y1+1 then
	begin vrastplot(S_dir,x0,y1,0); go to finished;
	end
else y0←y0-1;
good_case: if y0≤y1 then vrastplot(N_dir,x0,y0,y1-y0)
	else vrastplot(S_dir,x0,y1,y0-y1);
finished: lastx←x0; lasty←y1;
end;

procedure hplot(integer x0,y0,x1) # current pen applied to $(x0:x1,y0)$;
begin comment This procedure uses the current pen or eraser to update the raster
in a horizontal line going left or right from $(x0,y0)$ to $(x0,y1)$ inclusive;
label plot_one # go here if $(x0,y0)$ must be fully plotted with interior points;
label good_case # go here if the same edge can be used for the whole plot;
label remove_one # go here to remove the first point $(x0,y0)$;
label finished # go here when the plotting is done;
integer p # pointer to proper edge of current pen;
if x0 min x1 < xrastmin or x0 max x1 > xrastmax or y0<yrastmin or y0>yrastmax
	then drawerror("Curve out of range ("&cvs(x0)&":"&cvs(x1)&","&cvs(y0)&")");
if plotted then
	begin integer i,j;
	i←x0-lastx; j←y0-lasty;
	if abs(i)>1 or abs(j)>1 then go to plot_one;
	case 3*i+j+4 of begin
	[0] begin vrastplot(SW_dir,x0,y0,0); go to remove_one end;
	[1] if x1<x0 then go to good_case else begin hrastplot(W_dir,x0,y0,0);
	go to remove_one end;
	[2] begin vrastplot(NW_dir,x0,y0,0); go to remove_one end;
	[3] begin vrastplot(S_dir,x0,y0,0); go to remove_one end;
	[4] go to good_case;
	[5] begin vrastplot(N_dir,x0,y0,0); go to remove_one end;
	[6] begin vrastplot(SE_dir,x0,y0,0); go to remove_one end;
	[7] if x1≥x0 then go to good_case else begin hrastplot(E_dir,x0,y0,0);
	go to remove_one end;
	[8] begin vrastplot(NE_dir,x0,y0,0); go to remove_one end;
	else confusion
	  end;
	end
else plotted←true;
plot_one: fullrastplot(x0,y0);
remove_one: if x0=x1 then go to finished;
if x0<x1 then x0←x0+1
else if x0=x1+1 then
	begin hrastplot(W_dir,x1,y0,0); go to finished;
	end
else x0←x0-1;
good_case: if x0≤x1 then hrastplot(E_dir,x0,y0,x1-x0)
	else hrastplot(W_dir,x1,y0,x0-x1);
finished: lastx←x1; lasty←y0;
end;
comment The routine that plots a cubic;

procedure spdraw(real rw,rx0,rx1,rx2,rx3,ry0,ry1,ry2,ry3);
begin comment This procedure plots the cubic curves defined by 
$x(t)= \\{rx0}+\\{rx1}\cdot t+\\{rx2}\cdot t↑2+\\{rx3}\cdot t↑3$ and
$y(t)= \\{ry0}+\\{ry1}\cdot t+\\{ry2}\cdot t↑2+\\{ry3}\cdot t↑3$ for $0≤t≤1$,
where the pen size is \\{rw}.

First the interval is subdivided, if necessary, so that the relevant cubics are
monotonic and so that we don't have both $|x↑\prime(t)| > |y↑\prime(t)|$ and
$|x↑\prime(t)| < |y↑\prime(t)|$ in the same interval. The variable \\{octant} will
be used to classify the subintervals according to the following scheme:

	octant=0	0 ≤ y'(t) ≤ x'(t)
	octant=1	0 ≤ y'(t) ≤ -x'(t)
	octant=2	0 ≤ -y'(t) ≤ x'(t)
	octant=3	0 ≤ -y'(t) ≤ -x'(t)
	octant=4	0 ≤ x'(t) ≤ y'(t)
	octant=5	0 ≤ -x'(t) ≤ y'(t)
	octant=6	0 ≤ x'(t) ≤ -y'(t)
	octant=7	0 ≤ -x'(t) ≤ -y'(t)

When $\\{octant}=0$, for example, the drawing routine makes use of the facts that
$x(t)$ and $y(t)$ are nondecreasing and that $y(t)$ is growing less rapidly than
$x(t)$. The subintervals are specified by an array of points
	1 = tt[0] > tt[1] > ... > tt[m] = 0
where $\\{tt}[i-1]≥\\{tt}[i]+\\{eps}$ for $1≤i≤m$.

Since METAFONT works with a discrete raster, it gives careful attention to the set
of integer points $(k,l)$ that are plotted. The general idea is to plot $(k,l)$ if,
for some $0≤t≤1$, we have $k=\round\biglp x(t)\bigrp$ and 
$l=\round\biglp y(t)\bigrp$ and $|x(t)-k|+|y(t)-l|<{1\over2}$. In other words,
we plot $(k,l)$ when the curve hits the diamond-shaped region
		|x-k| + |y-l| < 1/2.
We also plot $(k,l)$ in certain boundary cases when the curve just happens to go
between diamonds, or when it appears to do so as a consequence of the limited
accuracy of floating-point arithmetic. A fairly arbitrary decision is made in
such boundary cases, but this procedure preserves important symmetry: If $x(t)$
were to be replaced by $n-x(t)$ for some integer $n$, the set of plotted points
$(k,l)$ would be replaced by the set of points $(n-k,l)$, except in the case that
$x(t)$ is constant and equal to an integer plus 1/2. (In the latter case such
symmetry is obviously impossible.) Similarly there would be reflective symmetry
if $y(t)$ were changed to $n-y(t)$, or if $x(t)$ and $y(t)$ were interchanged.
This symmetry is a simple consequence of the fact that octants are used to govern
the plotting: In each subinterval the curves are reduced to the case
$\\{octant}=0$ by performing the operations $x(t)←-x(t)$ and/or $y(t)←-y(t)$
and/or $x(t)↔y(t)$. There operations are applied later in reverse order when
drawing the actual curve in the actual octant.

A correction is added to $x(t)$ and/or $y(t)$ to compensate for possible
off-centeredness of a pen whose width and/or height is even. The routine is
designed to plot a straight horizontal or vertical sequence of points by making
a single call to \\{hplot} or \\{vplot}. Only the edges of the pen are plotted,
corresponding to increasing values of $t$.
;
integer octant # one of eight ways the curve may be varying (see above);
real a0,a1,a2,a3,b0,b1,b2,b3 # coefficients of transformed $x(t)$ and $y(t)$;
integer m # the number of subintervals;
integer w # pen size;

procedure ttinsert(real t);
begin comment This procedure inserts $t$ into the \\{tt} table, unless $t$ is
out of bounds or within \\{eps} of one the entries already present;
integer i;
if t>1.0-eps or t<eps then return;
i←1; while t≤tt[i] do i←i+1;
if t>tt[i-1]-eps or t<tt[i]+eps then return;
while i≤m do
	begin t↔tt[i]; i←i+1;
	end;
m←i; tt[m]←0.0;
end;

procedure split(real a,b,c);
begin comment This procedure calls $\\{ttinsert}(t)$ for all $t$ such that the
polynomial $3at↑2+2bt+c$ changes sign at $t$;
a←3.0*a;
if abs(a)≥eps then
	begin real d,r;
	if (d←b↑2-a*c)>eps then
		begin if b≥0 then r←(-b-sqrt(d))/a
		else r←(-b+sqrt(d))/a;
		ttinsert(r); ttinsert((c/a)/r);
		end;
	end
else if abs(b)≥eps then ttinsert(-.5*c/b);
end;

recursive procedure initlr(integer l0,l1);
begin comment This subprocedure is used to fill the \\{xl}, \\{xr}, \\{tl}, and
\\{tr} arrays in positions between $l0$ and $l1$, according to the conventions
described below. It is assumed that $\\{xr}[l0]$, $\\{tr}[l0]$, $\\{xl}[l1]$, and
$\\{tl}[l1]$ are already known;
integer k,l; real t0,t1,t;
if l1≤l0+1 then return;
t0←tr[l0]; t1←tl[l1];
loop	begin t←.5*(t0+t1);
	if t-t0≤eps then drawerror("Curve too wild");
	k←((a3*t+a2)*t+a1)*t+a0;
	l←((b3*t+b2)*t+b1)*t+b0;
	if l=l0 then
		begin tr[l]←t; t0←t; xr[l]←k;
		end
	else if l=l1 then
		begin tl[l]←t; t1←t; xl[l]←k;
		end
	else done;
	end;
tl[l]←tr[l]←t;
xl[l]←xr[l]←k;
initlr(l0,l);
initlr(l,l1);
end;

comment The \\{spdraw} procedure begins here;
if plottrace then print(nextline);
tt[0]←1.0; tt[1]←0.0; m←1;
split(rx3,rx2,rx1);
split(ry3,ry2,ry1);
split(rx3+ry3,rx2+ry2,rx1+ry1);
split(rx3-ry3,rx2-ry2,rx1-ry1);
comment Now each subinterval belongs to a single octant;
w←rw+.5 # round the pen size;
setuppen(w);
define corr(val,size)=⊂if size land 1 then val else val-.5⊃;
case curpen of begin
[cpen] begin rx0←corr(rx0,w); ry0←corr(ry0,w) end;
[hpen] begin rx0←corr(rx0,w); ry0←corr(ry0,hpenht) end;
[vpen] begin rx0←corr(rx0,vpenwd); ry0←corr(ry0,w) end;
[lpen] ry0←corr(ry0,lpenht);
[rpen] ry0←corr(ry0,rpenht);
[spen] begin rx0←rx0-sxcorr; ry0←ry0-sycorr end;
[epen] begin rx0←rx0-excorr; ry0←ry0-eycorr end;
else confusion
  end;

while m>0 do
	begin comment Now we draw the curve between tt[m] and tt[m-1];
	real t0,t1,t,xprime,yprime,xi,eta;
	integer x0,x1,ax,bx,ay,by,l;
	t0←tt[m]; t1←tt[m-1]; m←m-1;
	t←.5*(t0+t1);
	xprime←(3*rx3*t+2*rx2)*t+rx1;
	if xprime≥0 then
		begin a0←rx0; a1←rx1; a2←rx2; a3←rx3; octant←0;
		end
	else	begin comment $x(t)$ decreasing, negate the sign;
		a0←-rx0; a1←-rx1; a2←-rx2; a3←-rx3; octant←1;
		end;
	yprime←(3*ry3*t+2*ry2)*t+ry1;
	if yprime≥0 then
		begin b0←ry0; b1←ry1; b2←ry2; b3←ry3;
		end
	else	begin comment $y(t)$ decreasing, negate the sign;
		b0←-ry0; b1←-ry1; b2←-ry2; b3←-ry3; octant←octant+2;
		end;
	if abs(yprime)>abs(xprime) then
		begin comment $y$ growing faster thn $x$, interchange them;
		a0↔b0; a1↔b1; a2↔b2; a3↔b3; octant←octant+4;
		end;
	b0←b0+.5 # reduce rounding to truncation;

	xi←((a3*t1+a2)*t1+a1)*t1+a0;
	eta←((b3*t1+b2)*t1+b1)*t1+b0;
	x1←bx←xi; by←eta;
	if abs(xi-bx-1.0)+abs(eta-by-.5)≤.5 then bx←bx+1;
	xi←((a3*t0+a2)*t0+a1)*t0+a0;
	eta←((b3*t0+b2)*t0+b1)*t0+b0;
	x0←ax←xi; ay←eta;
	if abs(xi-ax)+abs(eta-ay-.5)≤.5 then ax←ax-1;
	comment Taking due account of whether or not this curve segment touches
	the diamonds at its left and right ends, we will plot points $(k,l)$ for
	$\\{ax}<k≤\\{bx}$. The relevant $l$ values are $\\{ay}≤l≤\\{by}$. For
	each $l$ in this range we will maintain table entries $\\{tl}[l]$ and
	$\\{tr}[l]$, the least and greatest $t$ values known such that
	$\round\biglp y(t)\bigrp=l$. We also will store the corresponding
	truncated $x(t)$ values $\\{xl}[l]$ and $\\{xr}[l]$. Our goal is to
	refine our knowledge about the cubics so that
		xr[l] = xl[l+1]  for  ay ≤ l < by.
	Then we shall plot $(k,l)$ for
		xl[l] < k ≤ xr[l],  ay ≤ l ≤ by.
	The boundary values $\\{xl}[{ay}]$ and $\\{xr}[{by}]$ are set to
	$\\{ax}$ and $\\{bx}$, respectively.
	;
	if ay<pmin or by>pmax then drawerror("Curve out of range ("&cvs(ay)&
		":"&cvs(by)&")");
	xr[ay]←x0; tr[ay]←t0;
	xl[by]←x1; tl[by]←t1;
	xl[ay]←ax; xr[by]←bx # note order of assignment in case $\\{ay}=\\{by}$;
	initlr(ay,by);
	for l←ay thru by do
		begin if l≠by and xr[l]<xl[l+1] then
			begin real t0,t1,t; integer k0,k1,k,ll;
			t0←tr[l]; t1←tl[l+1]; k0←xr[l]; k1←xl[l+1];
			loop	begin t←.5*(t0+t1);
				if t-t0≤eps then
					begin k1←k1-1;
					if k0=k1 then done;
					drawerror("Curve too wild");
					end
				else	begin k←((a3*t+a2)*t+a1)*t+a0;
					ll←((b3*t+b2)*t+b1)*t+b0;
					if ll=l then
						begin t0←t; k0←k;
						end
					else	begin t1←t; k1←k;
						end;
					if k0=k1 then done;
					end;
				end;
			xr[l]←k0; xl[l+1]←k1;
			comment There's no need to reset $\\{th}[l]$, $\\{tl}[l+1]$;
			end;
		if xl[l]<xr[l] then
			begin comment plotting $(k,l)$ for all appropriate $k$;
			case octant of begin
			[0] hplot(xl[l]+1,l,xr[l]);
			[1] hplot(-xl[l]-1,l,-xr[l]);
			[2] hplot(xl[l]+1,-l,xr[l]);
			[3] hplot(-xl[l]-1,-l,-xr[l]);
			[4] vplot(l,xl[l]+1,xr[l]);
			[5] vplot(-l,xl[l]+1,xr[l]);
			[6] vplot(l,-xl[l]-1,-xr[l]);
			[7] vplot(-l,-xl[l]-1,-xr[l]);
			else confusion
			  end;
			end;
		end;
	end;
end;
procedure ccubics # compute the cubic splines;
begin comment This procedure uses the data in arrays \\{pointx}, \\{pointy},
\\{tanx}, and \\{tany} to compute the coefficients of cubic splines to be drawn;
integer j,k;
label trivial;
if npts=1 then go to trivial;

comment First we set up auxiliary difference tables;
for k←0 thru npts do
	begin xdel[k]←pointx[k+1]-pointx[k];
	ydel[k]←pointy[k+1]-pointy[k];
	del[k]←xdel[k]↑2+ydel[k]↑2;
	end;
comment Let $\Delta z = \Delta x + i\Delta y$ in the interval from the $k$th
point to the $(k+1)$st point. Then $\\{xdel}[k]=\Delta x$ and $\\{ydel}[k]=\Delta y$
and $\\{del}[k]=|\Delta z|↑2$. We will compute the tangent directions so that if the
cubic spline begins in the direction of $e↑{i\theta}\Delta z$ and ends in the
direction of $e↑{-i\varphi}\Delta z$ then $e↑{i\theta}=(\\{tanx}[k]+i\cdot
\\{tany}[k])|\Delta z|/\Delta z$ and $e↑{-i\varphi}=(\\{tanx}[k+1]+i\cdot
\\{tany}[k+1])|\Delta z|/\Delta z$;
for k←1 thru npts do
	begin real alpha,beta # desired tangent direction at current point;
	real gamma # $\sqrt{\\{alpha}↑2+\\{beta}↑2}$;
	if abs(tanx[k])>eps or abs(tany[k])>eps then
		begin alpha←tanx[k]; beta←tany[k] # direction specified by user;
		end
	else if del[k-1]≤eps then
		begin alpha←xdel[k]; beta←ydel[k] # first point of series;
		end
	else if del[k]≤eps then
		begin alpha←xdel[k-1]; beta←ydel[k-1] # last point of series;
		end
	else	begin comment Choose tangent direction determined by the circle
			passing through points $k-1$, $k$, and $k+1$;
		alpha←xdel[k-1]/del[k-1]+xdel[k]/del[k];
		beta←ydel[k-1]/del[k-1]+ydel[k]/del[k];
		end;
	gamma←sqrt(alpha↑2+beta↑2);
	if gamma≤eps then
		begin alpha←gamma←1.0; beta←0.0;
		end;
	tanx[k]←alpha/gamma; tany[k]←beta/gamma;
	comment Now $\\{tanx}[k]↑2+\\{tany}[k]↑2=1$;
	end;
for k←1 thru npts-1 do
	begin label linear # go here if the curve is linear not cubic;
	real tcos, tsin # $\\{tcos}+i\cdot\\{tsin}=e↑{i\theta}|\Delta z|$;
	real fcos, fsin # $\\{fcos}+i\cdot\\{fsin}=e↑{i\varphi}|\Delta z|$;
	real ctp # $\cos(\theta+\varphi)$;
	real ctp2 # $\cos\biglp(\theta+\varphi)/2\bigrp$;
	real delta # $\sqrt(\\{del}[k])$;
	real alpha # temporary storage;
	real r, s # "velocities";
	x[0,k]←pointx[k]; y[0,k]←pointy[k] # the constant coefficient;
	if del[k]≤eps then go to linear;
	tcos←xdel[k]*tanx[k]+ydel[k]*tany[k];
	tsin←-ydel[k]*tanx[k]+xdel[k]*tany[k];
	fcos←xdel[k]*tanx[k+1]+ydel[k]*tany[k+1];
	fsin←ydel[k]*tanx[k+1]-xdel[k]*tany[k+1];
	ctp←(tcos*fcos-tsin*fsin)/del[k];
	delta←sqrt(del[k]);
	if ctp<.99999 then
		begin ctp2←sqrt(.5*(0.0 max(ctp+1.0)));
		alpha←sqrt(8.0)/(sqrt(1.0-ctp)*(1.0+ctp2));
		r←abs(alpha*fsin); s←abs(alpha*tsin);
		end
	else	begin comment $\theta \approx -\varphi$, probably both zero;
		r←s←delta; comment changed from 2.0*delta (JDH, November 83);
		end;
	if r>maxvr*delta then
		begin if modtrace then print(nextline,
		"Spline velocity reduced between points ",
		indexname(pointi[k])," and ",indexname(pointi[k+1]),
		" (r =",cvf(r/delta),")");
		r←maxvr*delta;
		end;
	if s>maxvs*delta then
		begin if modtrace then print(nextline,
		"Spline velocity reduced between points ",
		indexname(pointi[k])," and ",indexname(pointi[k+1]),
		" (s =",cvf(s/delta),")");
		s←maxvs*delta;
		end;
	if r<minvr*delta then
		begin if modtrace then print(nextline,
		"Sharp turn suppressed between points ",
		indexname(pointi[k])," and ",indexname(pointi[k+1]),
		" (r =",cvf(r/delta),")");
		r←minvr*delta;
		end;
	if s<minvs*delta then
		begin if modtrace then print(nextline,
		"Sharp turn suppressed between points ",
		indexname(pointi[k]), " and ",
		indexname(pointi[k+1]),
		" (s =",cvf(s/delta),")");
		s←minvs*delta;
		end;
	comment the desired cubic for $x$ is now
		pointx[k] + t\cdot xdel[k] + t(1-t)((1-t)α-tβ)
	where $α$ and $β$ are the desired $x$-derivatives minus $\\{xdel}[k]$:
		α = r\cdot tanx[k] - xdel[k],  β = s\cdot tanx[k+1] - xdel[k];
	x[1,k]←r*tanx[k]; alpha←x[1,k]-xdel[k];
	x[3,k]←alpha+s*tanx[k+1]-xdel[k] # $α+β$;
	x[2,k]←-x[3,k]-alpha # $-2α-β$;
	y[1,k]←r*tany[k]; alpha←y[1,k]-ydel[k];
	y[3,k]←alpha+s*tany[k+1]-ydel[k] # $α+β$;
	y[2,k]←-y[3,k]-alpha # $-2α-β$;
	continue;
	linear: x[1,k]←xdel[k]; y[1,k]←ydel[k];
	x[2,k]←x[3,k]←y[2,k]←y[3,k]←0.0;
	end;

trivial: x[0,npts]←pointx[npts]; y[0,npts]←pointy[npts];
x[1,npts]←x[2,npts]←x[3,npts]←y[1,npts]←y[2,npts]←y[3,npts]←0.0;

comment Finally we rotate the results;
if xxtr≠1.0 or xytr or xtr or yxtr or yytr≠1.0 or ytr then
	begin for k←1 thru npts do for j←0 thru 3 do
		begin real temp;
		temp←xxtr*x[j,k]+xytr*y[j,k];
		y[j,k]←yxtr*x[j,k]+yytr*y[j,k];
		if j=0 then
			begin y[0,k]←y[0,k]+ytr; temp←temp+xtr;
			end;
		x[j,k]←temp;
		end;
	end;
end;
procedure filldraw(real pensize);
begin comment This procedure will use the spline curves defined in arrays
\\x and \\y, and the spline curves defined in arrays \\{xx} and \\{yy}, to
draw a series of curves that will fill in between these extremes;
integer w # rounded pen size;
integer k # runs through the specified pairs of points;
integer n # one less than the number of curves to be drawn;
integer j # runs from 0 to $n$;

w←pensize+.5;
setuppen(w); n←1 # always draw at least the two extreme curves;
for k←1 thru npts do
	begin comment Figure the number of points needed to go from $(x,y)$ to
	$(\\{xx},\\{yy}) in a straight line with the specified pen, at point $k$;
	integer m;
	real alpha,beta;
	real d # maximum distance we could theoretically put between overlapping
		pens in the desired direction, if resolution were infinite;
	alpha←abs(x[0,k]-xx[0,k]); beta←abs(y[0,k]-yy[0,k]);
	if alpha>eps or beta>eps then
		begin case curpen of begin
		[cpen][spen][epen] d←sqrt(alpha↑2+beta↑2)/w;
		[hpen] d←sqrt((alpha/w)↑2+(beta/hpenht)↑2);
		[vpen] d←sqrt((alpha/vpenwd)↑2+(beta/w)↑2);
		[lpen] d←(alpha/w)max(beta/lpenht);
		[rpen] d←(alpha/w)max(beta/rpenht);
		else drawerror("No pen defined")
		  end;
		comment Working with infinite resolution, we could safely use
		$\lceil d\rceil +1$ points to connect $x[0,k]$ to $\\{xx}[0,k]$;
		m←d*safetyfactor + 1 # default safetyfactor is 2.0;
		if m>n then n←m;
		end;
	end;

comment We will use $n+1$ equally spaced spline curves;
for j←0 thru n do
	begin real alpha; alpha←j/n;
	for k←1 thru npts-1 do
		spdraw(w,x[0,k]+alpha*(xx[0,k]-x[0,k]),
				x[1,k]+alpha*(xx[1,k]-x[1,k]),
				x[2,k]+alpha*(xx[2,k]-x[2,k]),
				x[3,k]+alpha*(xx[3,k]-x[3,k]),
				y[0,k]+alpha*(yy[0,k]-y[0,k]),
				y[1,k]+alpha*(yy[1,k]-y[1,k]),
				y[2,k]+alpha*(yy[2,k]-y[2,k]),
				y[3,k]+alpha*(yy[3,k]-y[3,k]));
	end;

comment Finally we fill in the ends;
spdraw(w,x[0,1],xx[0,1]-x[0,1],0.0,0.0,y[0,1],yy[0,1]-y[0,1],0.0,0.0);
spdraw(w,x[0,npts],xx[0,npts]-x[0,npts],0.0,0.0,
	y[0,npts],yy[0,npts]-y[0,npts],0.0,0.0);
end;
comment The \\{drawit} procedure begins here;

integer k # runs through the specified points;
integer j # runs from 0 to 3 (for coefficients);

if drawtrace then
	begin integer wd,dg # width and number of digits printed;
	getformat(wd,dg);
	print(nextline,if ddrawit then "ddraw:  (" else "draw:   (",
		strpen,if eraser then "#)  " else ")   ");
	if not ddrawit then
		begin for k←1 thru npts-1 do
			print(if pointstab[k] then "#" else " ","   ...   ");
		if pointstab[npts] then print("#");
		print("     ");
		end;
	if ddrawit then
		begin print(nextline,"       ");
		for k←0 thru npts+1 do print(if dpnti[k]<0 then "        "
			else (indexname(dpnti[k])&"          ")[1 to 10]);
		setformat(10,4);
		print(nextline," x =");
		for k←0 thru npts+1 do print(dpntx[k]);
		print(nextline," y =");
		for k←0 thru npts+1 do print(dpnty[k]);
		print(nextline,"        tanx =");
		for k←1 thru npts do print(dtanx[k]);
		print(nextline,"        tany =");
		for k←1 thru npts do print(dtany[k]);
		setformat(wd,dg);
		end;
	print(nextline,"       ");
	for k←0 thru npts+1 do print(if pointi[k]<0 then "          "
		else (indexname(pointi[k])&"          ")[1 to 10]);
	setformat(10,4);
	if not ddrawit then
		begin print(nextline," w =");
		for k←0 thru npts+1 do print(pointw[k]);
		end;
	print(nextline," x =");
	for k←0 thru npts+1 do print(pointx[k]);
	print(nextline," y =");
	for k←0 thru npts+1 do print(pointy[k]);
	print(nextline,"        tanx =");
	for k←1 thru npts do print(tanx[k]);
	print(nextline,"        tany =");
	for k←1 thru npts do print(tany[k]);
	setformat(wd,dg);
	end;

plotted←false;
ccubics;
if ddrawit then
	begin 	pointw[1]←cursize max 1.0;
	for k←0 thru npts+1 do
		begin if k>0 then for j←0 thru 3 do
			begin xx[j,k]←x[j,k]; yy[j,k]←y[j,k];
			end;
		pointx[k]←dpntx[k]; pointy[k]←dpnty[k];
		pointi[k]←dpnti[k];
		tanx[k]←dtanx[k]; tany[k]←dtany[k];
		end;
	ccubics;
	filldraw(pointw[1]);
	end
else	begin real maxw,minw # maximum and minimum size of pen;
	integer mxw,mnw # rounded equivalents of \\{maxw} and \\{minw};
	maxw←minw←pointw[1];
	for k←2 thru npts do
		begin if pointw[k]<minw then minw←pointw[k]
		else if pointw[k]>maxw then maxw←pointw[k];
		end;
	mxw←maxw+.5; mnw←minw+.5;
	if mxw=mnw then
		begin comment The nice case (constant pen size);
		for k←1 thru npts-1 do
			spdraw(maxw,x[0,k],x[1,k],x[2,k],x[3,k],
				y[0,k],y[1,k],y[2,k],y[3,k]);
		end
	else	begin comment Varying pen size, reduce to ddraw;
		for k←0 thru npts do wdel[k]←pointw[k+1]-pointw[k];
		for k←1 thru npts do
			if pointstab[k] then tanw[k]←0.0
			else if del[k-1]≤eps then tanw[k]←wdel[k]
			else if del[k]≤eps then tanw[k]←wdel[k-1]
			else tanw[k]←(wdel[k-1]/del[k-1]+wdel[k]/del[k])/
				(1.0/del[k-1]+1.0/del[k]);
		for k←1 thru npts-1 do
			begin real alpha; w[1,k]←tanw[k];
			alpha←tanw[k]-wdel[k];
			w[3,k]←alpha+tanw[k+1]-wdel[k];
			w[2,k]←-w[3,k]-alpha;
			w[0,k]←pointw[k]-minw;
			end;
		w[0,npts]←pointw[npts]-minw;
		w[1,npts]←w[2,npts]←w[3,npts]←0.0;
		for k←1 thru npts do for j←0 thru 3 do
			begin case curpen of begin
			[hpen] begin xx[j,k]←x[j,k]-.5*w[j,k];
			x[j,k]←x[j,k]+.5*w[j,k];
			yy[j,k]←y[j,k] end;
			[vpen] begin xx[j,k]←x[j,k];
			yy[j,k]←y[j,k]-.5*w[j,k];
			y[j,k]←y[j,k]+.5*w[j,k] end;
			[lpen] begin xx[j,k]←x[j,k]-w[j,k];
			yy[j,k]←y[j,k] end;
			[rpen] begin xx[j,k]←x[j,k]+w[j,k];
			yy[j,k]←y[j,k] end;
			else drawerror("You can't vary the pen size with "
				&strpen)
			  end;
			end;
		filldraw(minw);
		end;
	end;

if not plotted then
	begin comment Single point or a curve that was too tiny;
	integer x0,y0,w;
	define ccor(val,size)=⊂if size land 1 then val+.5 else val⊃;
	w←pointw[1]+.5;
	setuppen(w);
	case curpen of begin
	[cpen] begin x0←ccor(x[0,1],w); y0←ccor(y[0,1],w) end;
	[hpen] begin x0←ccor(x[0,1],w); y0←ccor(y[0,1],hpenht) end;
	[vpen] begin x0←ccor(x[0,1],vpenwd); y0←ccor(y[0,1],w) end;
	[lpen] begin x0←x[0,1]+.5; y0←ccor(y[0,1],lpenht) end;
	[rpen] begin x0←x[0,1]+.5; y0←ccor(y[0,1],rpenht) end;
	[spen] begin x0←x[0,1]+.5-sxcorr; y0←y[0,1]+.5-sycorr end;
	[epen] begin x0←x[0,1]+.5-excorr; y0←y[0,1]+.5-eycorr end;
	else confusion
	  end;
	if x0<xrastmin or x0>xrastmax or y0<yrastmin or y0>yrastmax then
		drawerror("Point out of range");
	fullrastplot(x0,y0);
	end;

end_of_drawit: if drawdisplay then ddoutrast;
end;
comment the (temporary) routine for binary input;

internal integer bbuf # 0 or the number of rows in the next binput entry;

internal procedure binin # binary input to the raster;
begin integer y, yb, nn, xl, xr, rbits;
nn←wordin(bchan);
if ytr>yhigh then yhigh←ytr;
yb←ytr-bbuf+1;
if yb<ylow then ylow←yb;
if xtr<xrastmin+xpenmin or xtr+bitsperwd*nn>xrastmax+xpenmax or
	ytr>yrastmax+ypenmax or yb<yrastmin+ypenmin then
	begin error("binput out of range"); return;
	end;
xl←rcol(xtr); xr←rcol(bitsperwd*nn+xtr);
if xl<xleft then xleft←xl;
if xr>xright then xright←xr;
rbits←(xtr-xrastmin-xpenmin) mod bitsperwd;
for y←ytr step -1 until yb do
	begin integer i;
	for i←1 thru nn do
		begin integer w, x, xy;
		w←wordin(bchan);
		x←bitsperwd*(i-1)+xtr;
		xy←rloc(x,y);
		rast!gets!rast!lor!expr(xy,w lsh (-rbits))
		# rast[xy]←rast[xy] lor (w lsh (-rbits));
		rast!gets!rast!lor!expr(xy+rspan,w lsh (bitsperwd-rbits))
		# rast[xy+rspan]←rast[xy+rspan] lor (w lsh (bitsperwd-rbits));
		end;
	end;
bbuf←wordin(bchan) # get the next word (to see if it is the end of file);
if bbuf=0 then release(bchan);
end;
comment Now comes very system-dependent code for displaying on the user's screen;

IFDATADISC require "ddhdr.sai[hdr,he]" source_file; comment Hans Moravec's Datadisc package;
comment For documentation of Hans's routines see file \.{guide.gra[gra,hpm]};
external integer sline # location of row pointers in Hans's code;
saf integer array bp[0:7] # byte pointers for 4-bit groups;
integer bwd # buffer word used when converting from 36-bit to 32-bit format;
integer ddchan # Datadisc channel for output;
integer get,put # byte pointers for getting 4 bits and putting out 32;
preload_with true; safe boolean array ddnotready[0:0] # initialization needed;
internaldef ddxmin=-89,ddxmax=414,ddymin=-99,ddymax=380 # Datadisc window;
comment $\\{ddxmin}-1$ and \\{ddxmax} should be congruent to 18, modulo 36;
comment we must have xrastmin+xpenmin≤ddxmin, ddxmax≤xrastmax+xpenmax,
	yrastmin+ypenmin≤ddymin, ddymax≤yrastmax+ypenmax,
	ddxmax-ddxmin<504, ddymax-ddymin<480;
internaldef ddn=5 # printing is confined to this many lines at bottom of screen;

procedure initdd # do this before messing around with Datadisc routines;
begin integer k;
ddinit # Get Hans started;
screen(0,0,511,479);
ddchan←-1 # output goes to user's screen;
DEBUGONLY ddchan←gddchn(-1);
DEBUGONLY print(nextline,"Channel number ",cvos(ddchan));
DEBUGONLY erase(ddchan);
DEBUGONLY showa(ddchan);
for k←0 thru 7 do bp[k]←point(4,bwd,4*k+7);
DEBUGONLY if false then begin # omit the following line when debugging;
pppos(0,ddn*12-1) # set "page printer position";
DEBUGONLY end # omit the previous line when debugging;
end;

procedure cleardd # clear Hans's buffer;
begin if ddnotready[0] then return;
drken; rectan(0,0,511,479); liten;
end;

internal procedure ddoutrast # displays the raster on Datadisc screen;
begin comment The following code (suggested by Mike Farmwald) assumes that
$\\{bitsperwd}=36$ and refers to internal buffers of Hans's routines;
integer j,y,yl,yh,xl,xr,bytes;
	begin "am I a DD?"
	integer i;
	start_code  setom i; ttcall 6,i; end;  comment The GETLIN UUO;
	if i=-1 then i←0 # i=-1 means a detached job;
	if (i land '20000000000)=0 then
		begin error("Whoops, you need a Datadisc for display modes");
		control←control land lnot displaymodes;
		return;
		end;
	end "am I a DD?";
if ddnotready[0] then
	begin ddnotready[0]←false; initdd; cleardd;
	end;
xl←xleft max rcol(ddxmin); xr←xright min rcol(ddxmax);
yl←ylow max ddymin; yh←yhigh min ddymax;
bytes←9*(xr-xl+1) # number of 4-bit bytes to transmit;
for y←yl thru yh do
	begin integer xw; xw←xl*rspan+y;
	get←point(4,rast[xw],-1);
	put←point(32,memory[memory[location(sline)+ddymax-y]+location(dbuf)],-1);
		# the location of 32-bit pattern in Hans's buffer;
	bwd←0;
 	for j←0 thru bytes-1 do
		begin dpb(ildb(get),bp[j land 7]) # deposit 4 bits into \\{bwd};
		if (j land 7) = 7 then
			begin idpb(bwd,put); bwd←0;
			end;
		if (j mod 9) = 8 then
			begin xw←xw+rspan; 
			get←point(4,rast[xw],-1);
			end;
		end;
	if (bytes land 7) ≠ 0 then idpb(bwd,put) # deposit remaining bits;
	end;
dpyup(ddchan) # Hans's buffer to screen;
end; ENDDATADISC

IFADIS
IFTENEX
require "adis.sai" source_file;
ELSEC
require nextline&"Error: ADIS only implemented under TENEX!"&nextline message;
ENDC
ENDADIS
IFNOONLINE
procedure initdd;; procedure cleardd;; internal procedure ddoutrast;;
ENDNOONLINE
end